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Abstract 

The Mellin transform is usually applied in probability theory to the prod- 
uct of independent random variables. In recent times the machinery of the 
Mellin transform has been adopted to describe the Levy stable distributions, 
and more generally the probability distributions governed by generalized dif- 
fusion equations of fractional order in space and/or in time. In these cases 
the related stochastic processes are self-similar and are simply referred to as 
fractional diffusion processes. In this note, by using the convolution proper- 
ties of the Mellin transform, we provide some (interesting) integral formulas 
involving the distributions of these processes that can be interpreted in terms 
of subordination laws. 
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1. Introduction 



The role of the Mellin transform in probability theory is mainly related to 
the product of independent random variables: in fact it is well-known that 
the probability density of the product of two independent random variables 
is given by the Mellin convolution of the two corresponding densities. Less 
known is their role with respect to the class of the Levy stable distributions, 
that was formerly outlined by Zolotarev [3H] and Schneider [3D , see also [SB] . 
A general class of probability distributions (evolving in time), that includes 
the Levy strictly stable distributions, is obtained by solving, through the 
machinery of the Mellin transform, generalized diffusion equations of frac- 
tional order in space and/or in time, see [18\ 126] . In these cases the related 
stochastic processes turn out as self-similar and are referred to as fractional 
diffusion processes. 

In this note, after the essential notions and notations concerning the 
Mellin transform, we first show the role of the Mellin convolution between 
probability densities to establish subordination laws related to self-similar 
stochastic processes. Then, for the fractional diffusion processes we establish 
a sort of Mellin convolutions between the related probability densities, that 
can be interpreted as subordination laws. This is carried out starting for the 
representations through Merlin-Barnes integrals of the probability densities. 

We point out that our results, being based on simple manipulations, can 
be understood by non-specialists of transform methods and special functions; 
however they could be derived through a more general analysis involving the 
class of higher transcendental functions of Fox H type to which the prob- 
ability densities arising as fundamental solutions of the fractional diffusion 
equation belong. 

2. The Mellin transform 

The Mellin transform of a sufficiently well-behaved function f(x), x £ M + , 
is defined by 



when the integral converges. Here we assume f(x) € Li oc {]R + ) according to 
the most usual approach suitable for applied scientists. The basic proper- 
ties of the Mellin transform follow immediately from those of the bilateral 
Laplace transform since the two transforms are intimately connected. 
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(2.1) 
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Recently the theory of the Mellin transform has been the object 
of intensive researches by Professor Butzer and his associates, see e.g. 
[6j El [91 [lOl [11] ; in particular Butzer and Jansche [H [7] have introduced a 
theory independent from Laplace or Fourier transforms. 

The integral (2.1) defines the Mellin transform in a vertical strip in the 
s-plane whose boundaries are determined by the analytic structure of f(x) 
as x — > + and x — ► +oo. 

If we suppose that 

(0(i^- £ ) as x^0+, 

m = (2.2) 

I O (x 72 e ) as x — > +oo , 

for every (small) e > and 71 < 72, the the integral (2.1) converges abso- 
lutely and defines an analytic function in the strip 71 < 5Rs < 72 . This strip 
is known as the strip of analyticity of M {f(x); s} = f*(s) 

The inversion formula for (2.1) follows directly from the corresponding 
inversion formula for the bilateral Laplace transform. We have 

M- 1 {r(s);x} = f(x) = — / r(s)x- s ds, 7i<7<72, (2.3) 

at all points x > where f(x) is continuous. 

Let us now consider the most relevant Operational Rules. Denoting by 
the juxtaposition of a function f(x) on x > with its Mellin transform 

f*(s) , we have 

x a f(x) f*(s + a), a eW, 
f(x b ) ^ ±f*(s/b), beW, 6/0, 

f(cx) J^U c' s f*(s), ceM, c>0, 

from which 

x a f(cx b ) A ^" (S+a)/ "r(^) • (2-4) 
Furthermore we have 

h(x)= / o °°/(|) 5(0 f ^ /W« = *», (2-5) 
which is known as the Mellin convolution formula. 
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3. Subordination in stochastic processes via Mellin 
convolution 

In recent years a number of papers have appeared where explicitly or im- 
plicitly subordinated stochastic processes have been treated in view of their 
relevance in physical and financial applications, see e.g. [TJ [2J [3j EQl [3T1 (Ml 
[35| I3"6j [37] and references therein. Historically, the notion of subordination 
was originated by Bochner, see [U [5]. In brief, according to Feller [T7], a 
subordinated process X(t) = Y(T(t)) is obtained by randomizing the time 
clock of a stochastic process Y{t) using a new clock T(t), where T(t) is a 
random process with non-negative independent increments. The resulting 
process Y(T(t)) is said to be subordinated to Y(t), called the parent process, 
and is directed by T(t) called the directing process. The directing process 
is often referred to as the operational timqj. In particular, assuming Y(t) 
to be a Markov process with a spatial probability density function (pdf) of 
x, evolving in time r, q T (x) = q(x;r), and T(r) to be a process with non- 
negative independent increments with pdf of r depending on a parameter t, 
ut(r) = u(r;t), then the subordinated process X{t) = Y(T{t)) is governed 
by the spatial pdf of x evolving with t, pt(x) = p{x\ t), given by the integral 
representation 

poo 

Pt(x) = q T (x) u t (r) dr . (3.1) 
J o 

If the parent process Y(t) is self-similar of the kind that its pdf q T (x) is 
such that 

q T (x) = q(x; r) = q (jA , 7 > , (3.2) 
then Eq. (3.1) reads, 

Herewith we show how to interpret Eq. (3.3) in terms of a special convolution 
integral in the framework of the theory of the Mellin transform. Later, in 
the next sections, we shall show how to use the tools of the Mellin-Barnes 
integral and Mellin transform to treat the subordination for the class of 

1 ADDED NOTE (2007). In a recent paper by R. Gorenflo, F. Mainardi and A. 
Vivoli: Continuous time random walk and parametric subordination in fractional dif- 
fusion (http://arXiv.org/abs/cond-mat/0701126), to be published in Chaos, Solitons 
& Fractals (2007), the Authors have given an interesting alternative interpretation of 
subordination, that they call parametric subordination. 
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self-similar stochastic process, which are governed by fractional diffusion 
equations. 

Let X\ and X2 be two real independent random variables with pdfs 
Pi{x\) and P2{x2) respectively, with x\ G M and X2 £ Mq. The joint 
probability is 

P*(,x 1 ,x 2 ) = pi{x 1 )p 2 {x 2 ) ■ (3.4) 

Denoting by X the random variable obtained by the product of X\ and X.J , 
i.e. x = x\ x\ , and carrying out the transformation 

f xi = x/t~< , 

[X2 = T, 

we get the identity 



r) dx tir = pi(x/r 7 ) ^2(1") J dxdr , 



(3.5) 
(3.6) 



where 



J 





dx\ 




dx 


dx 2 




dx 2 


8x2 




dx 


8x2 









JX 
T 7+l 
1 



is the Jacobian of the transformation (3.5). Noting that J 
integrating (3.6) in dr we finally get the pdf of X, 



(3.7) 



1/t 7 and 



p(x) 



Oil 



Pi 



^)P2(.t)-, xeM. 



(3.8) 



For 7 = 1, by comparing with Eq. (2.5), we recover the well known result 
that the probability density of the product of two independent random vari- 
ables is given by the Mellin convolution of the two corresponding densities. 

We now adapt Eq. (3.8) to our subordination formula (3.3) by identifying 
Pi,P2 with q T and ut, respectively. We can now interpret the subordination 
formula (3.3) as follows. The pdf of the subordinated process X, pt(x), turns 
out to be the pdf of the product of the independent random variables X q 
and Xu 1 distributed according to q T {x q ) and ut(x u ), respectively. 



4. Fractional diffusion equations and probability 
distributions 

An interesting way to generalize the classical diffusion equation 
d 2 d 

— — Ku(x,t) = —u(x,t), — 00 < 2; < +00 , i>0, (4.1) 
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is to replace in (4.1) the partial derivatives in space and time by suitable 
linear integro-differential operators, to be intended as derivatives of non 
integer order, that allow the corresponding Green function (see below) to 
be still interpreted as a spatial probability density evolving in time with an 
appropriate similarity law. 

The Space- Time Fractional Diffusion Equation 

Recalling the approach by Mainardi, Luchko and Pagnini in [26], to 
which we refer the interested reader for details, it turns out that this general- 
ized diffusion equation, that we call space-time fractional diffusion equation, 
is 

x Dq u(x,t) = t D* u(x,t) , -oo<x<+oo, t>0, (4.2) 
where the a , 6 , (3 are real parameters restricted as follows 

0<a<2, \e\ < min{a,2 - a} , (4.3) 

< f3 < I or 1< p < a < 2 . (4.4) 

Here x Dq and tD* are integro-differential operators, the Riesz-Feller space 
fractional derivative of order a and asymmetry 6 and the Caputo time frac- 
tional derivative of order /3, respectively. 

The relevant cases of the space-time fractional diffusion equation (4.2) 
include, in addition to the standard case of normal diffusion {a = 2, (3 = 1}, 
the limiting case of the D'Alembert wave equation {a = 2,(3 = 2}, the 
space fractional diffusion {0 < a < 2, [3 = 1}, the time fractional diffusion 
{a = 2,0 < (3 < 2}, and the neutral fractional diffusion {0 < a = (3 < 2}. 
When 1 < (3 < 2 we speak more properly about the fractional diffusion-wave 
equation in that the corresponding equation governs intermediate phenom- 
ena between diffusion and wave processes. 

Let us now resume the essential definitions of the fractional derivatives 
in (4.2)-(4.4) based on their Fourier and Laplace representations. 

By denoting the Fourier transform of a sufficiently well-behaved (gener- 
alized) function f(x) as /(k) = J^{/(x); k} = /+~ e +iKX f(x) dx , k € M , 
the Riesz-Feller space fractional derivative of order a and skewness 8 is 
defined as 

^{*0?/(*);k} = -V£(«)/(«)> 

(4.5) 

V£( K ) = | K |<V(signK)07r/2 5 <a<2, \9\ < min {a, 2 - a} . 
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In other words the symbol of the pseudo-differential operator x Dg is the 
logarithm of the characteristic function of the generic strictly stable proba- 
bility density according to the Feller parameterization |16[ [T7] , as revisited 
by Gorenflo and Mainardi [21]. For this density we write 



T 



exp 



(4.6) 



where a is just the stability exponent (0 < a < 2) and is a real parameter 
related to the asymmetry (\0\ < min{a,2 — a}) improperly called skewness. 

By denoting the Laplace transform of a sufficiently well-behaved (gen- 
eralized) function f(t) as f(s) = C {f(t); s} = J °° e ~ st f(t) alt , R (s) > a f , 
the Caputo time fractional derivative of order (3 (m — 1 < (3 < m , mG M) 
is defined through^! 



m— 1 



k=0 



This leads to define, see e.g. [20|. [32] . 



1 



* /( m )(r)dr 



r(m-P)J (t 
d m 



)/3+l—m ' 



m 



1< /3 < m. (4.7) 



1 < P <m, 



(4.8) 



T/te Green Function 

When the diffusion equations (4.1), (4.2) are equipped by the initial and 
boundary conditions 

ti(x,0 + ) = (p(x) , u(±oo,t)=0, (4.9) 

their solution reads u(x, t) = G(£, t) ip(x — d£ , where G(x, t) denotes 
the fundamental solution (known as the Green function) corresponding to 
(p(x) = 5(x), the Dirac generalized function. We note that when 1 < (3 < 2 

2 The reader should observe that the Caputo fractional derivative introduced in 
[121 1131 114] represents a sort of regularization in the time origin for the classical Riemann- 
Liouville fractional derivative see e.g. [201132] . We note that the Caputo fractional deriva- 
tive coincides with that introduced (independently of Caputo) by Djrbashian & Nersesian 
[15] . which has been adopted by Kochubei, see e.g. [231 124] for treating initial value prob- 
lems in the presence of fractional derivatives. In [11] the authors have pointed out that 
such derivative was also considered by Liouville himself, but it should be noted that it 
was disregarded by Liouville who did not recognize its role. 
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we must add a second initial condition ut(x,0 + ) = ip(x), which implies 
two Green functions corresponding to {u(x,0 + ) = 5(x), ut(x,0 + ) = 0} and 
{u(x,0 + ) = 0, u t (x,0 + ) = S(x)}. Here we restrict ourselves to consider the 
first Green function because only for this it is legitimate to demand it to be 
a spatial probability density evolving in time, see below. 

It is straightforward to derive from (4.2) the composite Fourier-Laplace 
transform of the Green function by taking into account the Fourier transform 
for the Riesz-Feller space fractional derivative, see (4.5),, and the Laplace 
transform for the Caputo time fractional derivative, see (4.7). We have, see 
[26] 

& M ~ir£m- (4 - 10) 

By using the known scaling rules for the Fourier and Laplace transforms, 
we infer without inverting the two transforms, 

G%{x,t)=t-~<K%{x/t~<), 7 = /?/«, (4.11) 

where the one- variable function K e a * is the reduced Green junction and x/t 1 
is the similarity variable. We note from 

G^(0, S ) = l/ S ^G^(0,t) = l, 
the normalization property 

30 r+oo 

G a p{x,t)dx = / K a p(x) dx = 1 , 

<J — oo J— oo 

and, from 



<( K )=<(-k)=^(-k), 

the symmetry relation 

K%{-x)=K-°p{x), 

allowing us to restrict our attention to x > . 

For 1 < (5 < 2 we can show, see e.g. [28], that the second Green function 
is a primitive (with respect to the variable t) of the first Green function 
(4.11), so that, being no longer normalized in M, cannot be interpreted as 
a spatial probability density. 

When a = 2 (9 = 0) and (3 = 1 the inversion of the Fourier-Laplace 
transform in (4.10) is trivial: we recover the Gaussian density, evolving in 
time with variance a 2 = 2t, typical of the normal diffusion, 

G% fl (x,t) = ^^exp (-x 2 /(At)) , xGR, t>0, (4.12) 
which exhibits the similarity law (4.11) with 7 = 1/2. 
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For the analytical and computational determination of the reduced Green 
function we restrict our attention to x > because of the symmetry relation. 
In this range Mainardi, Luchko and Pagnini |26j have provided the Mellin- 
Barnea^l integral representation 

& (x)- 1 1 r io ° r(f)r(i-*)r(i- a ) x s ds 

a >P K ' ax2mj 1 - ioo r(l-|s)r(ps)r(l-ps) 2a ' 

(4.13) 

where 7 is a suitable real constant. 

The Space Fractional Diffusion : {0 < a < 2 , (3 = 1} 

In this case we recover the class L a (x) of the strictly stable (non- 
Gaussian) densities exhibiting fat tails (with the algebraic decay propor- 
tional to |x|~(° +1 )) and infinite variance, 

1 \ r f \ 11 P +lQ ° r(s/a) T(l - s) Sj a- 6 

KL(x) = L e a (x) = — -— / 7^ x ds ^ P ~ 



ax2iriJ 7 -i QO T(ps)T(l — ps) ' 2a ' 

(4.14) 

where < 7 < min{a, 1} . 

A stable pdf with extremal value of the skewness parameter is called 
extremal. One can prove that all the extremal stable pelfs' with < a < 1 
are one-sided, the support being Mq if 6 = —a , and Mq if 9 = +a . The 
one-sided stable pdf's with support in Mq can be better characterized by 
their (spatial) Laplace transform, which turn out to be 



oo 



Q: 



La a {s):= e~ bX L- a {x)dx = e~ b , (s) > , < a < 1 . (4.15) 







In terms of Mellin-Barnes integral representation we have 
1 1 f +>» r( s /a) ^ s 



La a (x) = — — -^V^ds, 0< 7 <a<l. (4.16) 

ax 2m J^-ioo 1 (s) 

The Time Fractional Diffusion : {a = 2 , < (3 < 2} 

In this case we recover the class of the Wright type densities exhibiting 
stretched exponential tails and finite variance proportional to t@ , 

1 1 1 n+ioo TVl _ <A 



3 The names refer to the two authors, who in the first 1910's developed the theory of 
these integrals using them for a complete integration of the hypergeometric differential 
equation. However, these integrals were first used by S. Pincherle in 1888. For a revisited 
analysis of the pioneering work by Pincherle (1853-1936, Professor of Mathematics at the 
University of Bologna from 1880 to 1928) we refer to Mainardi and Pagnini [27] . 
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where < 7 < 1 . As a matter of fact ^Mpi^ix) turns out to be a symmetric 
probability density related to the transcendental function M u (z) defined for 
any v £ (0, 1) and Vz G (E as 

00 (- 7 \ n 1 00 C—^n-l 

= E , rr rn — vi = - E rw r M sin (^)- ( 4 - 18 ) 

' n! 1 — un + (1 — i/) 7r , n — 1)! 

n=0 L v ' J n=l v ' 

We note that M v {z) is an entire function of order p = 1/(1 — v) , that turns 
out to be a special case of the Wright function^. Restricting our attention 
to the positive real axis (r > 0) we have: the Laplace transform 

£{M u (r);s} = E u {-s) , K(s)>0, < u < 1 , (4.19) 

where E u is the Mittag-Leffler function, and the asymptotic representation 

M v (r) ~ AqY v ~ x 1 2 exp (-Y) , x -» 00 , 

(4.20) 

a result formerly obtained by Wright himself, and, independently, by 
Mainardi and Tomirotti [29J by using the saddle point method. Because 
of the above exponential decay, any moment of order 5 > — 1 for M u (r) is 
finite and given as 

'r S M v (r)dr = ^ S+1) S > -1 . (4.21) 
V(u5 + 1) 

In particular we get the normalization property in M + , Jq°° M v (r) dr = 1. In 
view of Eqs (4.11) and (4.21) the moments (of even order) of the fundamental 
solution G® g(x, t) turn out to be, for n = 0, 1, 2, . . . and t > , 

> G%(x, t) dx = tf* [°°x 2n M m (x) dx = + t*\ (4.22) 
3 J i(pn + l) 



4 The Wright function is defined by the series representation, valid in the whole complex 
plane, 

n — 

Then, M v {z) := <l > _ J/ ,i_ J/ (— 2:) with < ^ < 1 . The function M„(z) provides a generaliza- 
tion of the Gaussian and of the Airy function in that 



M 



1/2 (z) = 4= ex P (- ^/4) , M 1/3 (z) = 3 2/3 Ai (z/3 1/3 ) 
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We agree to call M„(r) (y G (0,1), r G Mq) the M- Wright function of 
order u, understanding that its half represents the spatial pdf corresponding 
to the time fractional diffusion equation of order 2v. Relevant properties 
of this function, see e.g. [19, 25, 26], are concerning the limit expression 
for P = 1, i.e. M\{r) = 5{r — 1) and its relation with the extremal stable 
densities, i.e. 

^7^ 1/ (^j)=^T^(^) > 0<„<1, OO, r>0. (4.23) 

We note that, in both limiting cases of space fractional (a = 2) and time 
fractional (/? = 1) diffusion, we recover the Gaussian density of the normal 
diffusion, for which 



2x 2iri J^-ioo r(l — s/2) 
L° 2 (x)= 1 -M 1/2 ( X ) = ^ / =e- x2 ^. 



(4.24) 



T/te Neutral Fractional Diffusion : {0 < a = [3 < 2} 

In this case, surprisingly, the corresponding (reduced) Green function 
can be expressed (in explicit form) in terms of a (non-negative) simple ele- 
mentary function, that we denote by N^(x) , as it is shown in [26j : 



I i p+ioo r (jL)r(l-iL) 



AT^( X ) = / ^x s ds 

ax 2-ki J-y-ioc T{p s) F{1 — ps) 

1 x^sinff (a -6)] 
" 7T 1 + 2x a cos[f (a - 6»)] + x 2a ' 



(4.25) 



T/ie Fractional Diffusion Processes 

The self-similar stochastic processes generated by the above probability 
densities evolving in time can be considered as generalizations of the stan- 
dard diffusion processes and therefore distinguished from it with the label 
"fractional". When < (3 < 1 random walk models can be introduced 
to generalize the classical Brownian motion of the standard diffusion, as it 
was investigated in a number of papers of our group see e.g. \21\ \22\ . In 
the case of space fractional diffusion we obtain a special class of Marko- 
vian processes, called stable Levy motions, which exhibit infinite variance 
associated to the possibility of arbitrarily large jumps (Levy flights). In 
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the case of time fractional diffusion we obtain a class of stochastic pro- 
cesses which are non-Markovian and exhibit a variance consistent with slow 
anomalous diffusion. For the general genuine space-time fractional diffusion 
(0<a<2,0</3<l)we generate a class of densities (symmetric or 
non-symmetric according to 6 = or 6 ^ 0) which exhibit fat tails with an 
algebraic decay cx . Thus they belong to the domain of attraction of 

the Levy stable densities of index a and can be referred to as fractional sta- 
ble densities. The related stochastic processes possess the characteristics of 
the previous two classes; indeed, they are non-Markovian (being < (3 < 1) 
and exhibit infinite variance associated to the possibility of arbitrarily large 
jumps (being < a < 2). 



5. Subordination for space fractional diffusion pro- 
cesses 

In the book by Feller, see [T7], at p. 176 we read: Let X and Y be indepen- 
dent strictly stable variables, with characteristic exponent a and (3 respec- 
tively. Assume Y to be a positive variable (whence (3 < 1). The product 
X Y l l a has a stable distribution with exponent a /3. 

In other words, this statement means that any strictly stable process (of 
exponent 7 = a ■ /3) is subordinated to a parent strictly stable process (of 
exponent a) and directed by an extremal strictly stable process (of exponent 
< f3 < 1). Feller's proof is vague being, as a matter of fact, limited 
to symmetric subordinated and parent stable distributions. Furthermore, 
the proof, scattered in several sections, is essentially based on the use of 
Fourier and Laplace transforms. Here we would like to make more precise 
the previous statement by Feller by considering the possibility of asymmetry 
characterized by the index 9 as previously explained and making use of the 
Mellin machinery outlined in Sections 2 and 3, and of the Mellin-Barnes 
integral representation (4.14). So, in virtue of the Mellin inversion formula 
(2.3) we can write for the generic strictly stable pdf 

LS m 1 r(-^)rn + (.-1)] «-e 

Let us now consider the evolution in time according to the space fractional 
diffusion equation by writing 

4(z;t):=G4 fl (x,t)=*- 1/a 4(^) • (5-2) 
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Then we prove the following 

Theorem. Let L^ p (x;t), L%? q (x;i) and L d /3 (x;t) be strictly stable densities 
with exponents a p ,a g , {3 and asymmetry parameters 9 p ,9 q ,6/3, respectively, 
such that 

< otp < 2 , \6 P \ < min{a p , 2 — a p } , 
< a q < 2 , \9 q \ < mm{a q , 2 — a q } , 
0<P<1, Op = -0, 
then the following subordination formula holds true for < x < 00 , 

roo 

L% (x; t) = / L*» (x; r) l/(t; *) ^ , ^ » p = (3a q , 9 P = (39 q . (5.3) 

J 

Because of the scaling property (5.2) of the stable pd/'s, we can alternatively 
state 



00 



The proof of Eq. (5.4) is a (straightforward) consequence of the previous 
considerations. By recalling the Mellin pairs for the involved stable densities 
(that can be easily obtained from Eq. (5.1)-(5.2) by adopting the correct 
parameters) and the scaling properties of the Mellin transform, see (2.4), we 
have 

/ t \ \ a , / 1 \- s 1 r (-—) t\i + (s - 1) 

t -l/a pL 6 p ( t -l/a p ( 1 \ 1 \ a p ' 



^Kt^pJ ' Kt^pJ a p r[l + p p ( S -l)}T[-p p (s-l)} 

(5.5) 

and 



, a T O0, b-, M 1-42I r ( "bis 1 + }) r t 1 + ( S t a *)] 



After some algebra we recognize 

M{t-^L% ^y,s} = M{bx a cL^(cx b );s}M{L%(x);s} (5.7) 
provided that 

0p = -f3, a = a q -l, b = a q , c = t" 1//3 , (5.8) 

and 

a p = j3a q , 9 P = [39 q . (5.9) 
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Recalling the Mellin convolution formula (2.5) we obtain from Eqs. (5.7)- 
(5.9) the integral representation 

and, by replacing £ — > T l l aq , we finally get Eq. (5.4). 

Taking into account the relationships in Eq. (5.9) we can point out some 
interesting subordination laws. In particular we observe that any symmetric 
stable distribution with exponent a p = a S (0, 2) (6 p = 9 = 0) is subordi- 
nated to the Gaussian distribution (a g = 2, 6 q = 0), see L 2 {x) in (4.24), 
through an extremal stable density of exponent (5 = a/2, that is 

roc , 
L° a (x;t)=J o L° 2 (x;T)L a ^ 2 (T;t)dT, 0<a<2. (5.11) 

Furthermore, by recalling the generalized Cauchy density of skewness (\9\ < 
1), see e.g. Eq. (4.9) in [26], 

rBt ^ 1 COs(#7T/2) . n . 1 

LAX) = — ■: — - — j—tk = — — -—-TT , \V < 1 , -oo < X < +oo , 

iV y vr [x + sm(6»7r/2)] 2 + [cos(6»vr/2)] 2 1 1 

(5.12) 

we note that any stable distribution with exponent a p = a £ (0, 1) and 
skewness \9 p \ = \9\ < a is subordinated to the generalized Cauchy distribu- 
tion with skewness \6 q \ = \6\/a < 1 denoted by L^ a (x), see (5.12), through 
an extremal density of exponent (3 = a, that is 

L e a {x;t) = f°° L 9 1 /a (x;T)L- a {T-t)dT, 0<a<l. (5.13) 
J o 

6. Subordination for time fractional diffusion pro- 
cesses 

For the M-Wright function we deduce from (4.17) the Mellin transform pair: 

r[i + i/(s- 1)] 

Let us now consider the evolution in time according to the corresponding 
time fractional diffusion equation by writing 



M v (x-t):=2Gl 2v (x,t)=t-»M v [-) . (6.2) 
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Then we prove the following 

Theorem. Let M v {x;t), M v (x;t) and Mp{x;t) be M -Wright functions of 
orders u,r),(3 £ (0, 1) respectively, then the following subordination formula 
holds true for < x < oo, 

M u (x,t) = / M„(x; t ) Mq(t; t) dr , with u = rj (3 . (6.3) 
Jo 

Because of the scaling property (6.2) of the M- Wright functions, we can 
alternatively state 

The proof of Eq. (6.4) is a (straightforward) consequence of the previous 
considerations. After some algebra we recognize 

M{t~ v M v (^j ;s} = M{bcx a Mp(cx b );s} M{M v (x);s} (6.5) 

provided that 

a = - - 1 , b=-, c = t~ u ^ , (6.6) 

7] J] 

and 

u = /3rj. (6.7) 

Recalling the Mellin convolution formula (2.5) we obtain from Eqs. (6.5)- 
(6.7) the integral representation 

^(?)-jT;^^(f)^(^)f- <"> 

and, by replacing £ ^ t v , we finally get Eq. (6.4). 

We note that for rj = 1/2 the corresponding M- Wright function reduces 
to twice the Gaussian density according to (4.24), so the subordination for- 
mula (6.3) reads 

roo 

M m (x; t)=2 L° 2 (x; r) M p (r; t) dr . (6.9) 

J 
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7. Subordination for space-time fractional diffusion 



For the space-time fractional diffusion we shall prove two relevant subordi- 
nation laws. For this purpose we first consider the evolution in time writing 



K° 



t) := G%(x, t) = t~^ K% . (7.1) 

Then the required subordination laws read 

/•CO 

K 9 a A x > t )=J L e a {x-T)M p (r-t)dr, < < 1 . (7.2) 

and 

K%{x;t) = N e a (x;T)Mp /a (r;t)dT, </?/«< 1 • (7.3) 



In order to prove the above laws we start from two relevant results in 
see (6.1), that we report below 

[ a He" 1 Af„ (O] ^ a (x/0 f , < p < 1 , 

^£ (7 - 4) 
[ J Q M p/a {i) N e a (x/0 -j, < H I a < 1 . 

In [26] the above identities have enabled us to extend the probability interpre- 
tation of the Green functions to the ranges {0 < a < 2} n {0 < (3 < 1} and 
{1 < (3 < a < 2} . Indeed, the formulae in (7.4) show the non-negativity of 
the Green functions of the genuine space-time fractional diffusion equation 
based on the non-negativity of the Green functions of the particular cases 
of space, time and neutral fractional diffusion. We note that the formulae 
were derived by using the Mellin-Barnes representation of the corresponding 
Green functions, a method akin with that of the Mellin transform as we have 
previously seen. Let us now prove Eq. (7.2). Because of the scaling proper- 
ties of the involved functions this is equivalent to prove, with < /3 < 1, 

that can be easily achieved from the first equation of (7.4), by replacing x 
with xjt^l a and making the change of variable £ = r 1//a /t^/ a in the integral. 
Similarly, to prove Eq. (7.3) we can verify, with < — < 1, 

In this case it suffices to replace x with x/t^^ a and set £ = r/t^/ a in the 
second equation of (7.4). 
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It is worth to note that whereas for the particular cases of space and time 
fractional diffusion the corresponding subordination formulas (5.3) and (6.3) 
involve functions of the same class, here the subordination formulas (7.2) 
and (7.3) involve functions of different classes. However, in all the cases 
we point out that the pdf of the directing process is a M- Wright function, 
including (5.3) if we take into account (4.23). 

If we like to use the terminology G a g(x]t) for the Green functions, the 
two subordination laws (7.2)-(7.3) can be resumed as follows 

{poo 
2 J G e ajl (x;T)G° 2j2p (r;t)dT, < < 1 , 
2 ^ G e a Jx; r) G% 2/3/q (t; t)dr, 0< 0/a < 1 . 

From the first equation in (7.7) we note that the Green function for the 
space-time fractional diffusion equation of order {a, /?} , with < a < 2 and 
< /3 < 1 , can be expressed in terms of the Green function for the space 
fractional diffusion equation of order a and the Green function for the time 
fractional diffusion equation of order 2/3 . 



8. Concluding discussion 

There are several ways of deriving subordination formulas for fractional 
diffusion processes. A natural way is to to start from an approximating 
continuous time random walk model and carry out an appropriate passage 
to the limit. As we cannot give here a comprehensive survey of achievements 
of contributors to this subject, let us only cite [U|30j for a concise description. 
In contrast to this "stochastic" way it was our aim to show how classes of 
subordination formulas can be found by purely analytic methods, namely by 
exploiting solution formulas for the fundamental solutions of the essential 
types of fractional diffusion equations, formulas expressing these solutions 
in form of Mellin-Barnes integrals (see [26]). Using the machinery offered 
by Mellin transform theory these formulas can be rewritten as integrals of 
products that in turn allow a probabilistic interpretation as subordination 
formulas. So, we testify for the fact that the fascinating field of fractional 
diffusion processes is not only interesting from the view-points of probability 
theory and statistical physics, but also from that of pure analysis where it 
is a playground for lovers of special functions and integral transforms. In 
particular, we also demonstrate that the Mellin transform is an important 
and useful integral transform in its own right. 
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